% 高斯高通滤波
fxy=imread("img/lena.bmp");
subplot(2,2,1),imshow(fxy),xlabel("(a)原始图像");

% 傅里叶幅度谱的对数图像显示
fxy=double(fxy);
Fuv=fft2(fxy);%二维傅里叶变换
ShiftFuv=fftshift(Fuv);
AbsShiftFuv=abs(ShiftFuv);%傅里叶幅度谱
LogAbsShiftFuv=log(AbsShiftFuv);%傅里叶幅度谱可视化
subplot(2,2,2),imshow(LogAbsShiftFuv,[]),xlabel("(b)傅里叶幅度谱对数图像");



% 截至频率D0=10时的滤波效果
[N1,N2]=size(ShiftFuv);
D0=10;
%中心点
n1=fix(N1/2);
n2=fix(N2/2);
for i=1:N1
    for j=1:N2
        d=sqrt((i-n1)^2+(j-n2)^2);
        Huv=1-exp(-d^2/(2*D0^2));% 高斯高通公式
        G(i,j)=Huv*ShiftFuv(i,j);
    end
end
% 反变换
G=ifftshift(G);
g=ifft2(G);
g=uint8(real(g));
subplot(2,2,3),imshow(g),xlabel("(c)滤波图像(D0=10)");

% 截至频率D0=20时的滤波效果
[N1,N2]=size(ShiftFuv);
D0=20;
%中心点
n1=fix(N1/2);
n2=fix(N2/2);
for i=1:N1
    for j=1:N2
        d=sqrt((i-n1)^2+(j-n2)^2);
        Huv=1-exp(-d^2/(2*D0^2));
        G(i,j)=Huv*ShiftFuv(i,j);
    end
end
% 反变换
G=ifftshift(G);
g=ifft2(G);
g=uint8(real(g));
subplot(2,2,4),imshow(g),xlabel("(d)滤波图像(D0=20)");